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ABSTRACT 

Galaxy cluster surveys offer great promise for measuring cosmological parameters, but survey analy- 
sis methods have not been widely studied. Using methods developed decades ago for galaxy clustering 
studies, it is shown that nearly exact likelihood functions can be written down for galaxy cluster 
surveys. The sparse sampling of the density field by galaxy clusters allows simplifications that are 
not possible for galaxy surveys. An application to counts in cells is explicitly tested using cluster 
catalogs from numerical simulations and it is found that the calculated probability distributions are 
very accurate at masses above several times 10^'^h~^MQ at z = and lower masses at higher redshift. 

Subject headings: galaxy clusters; cosmology 



1. INTRODUCTION 

The idea of using cluster surveys as probes of 
dark energy has generated considerable interest of late 
llCarl strom et al. 2002; Haiman et al. 2001; Wellcr ct al. 
pool l20b.lfc iMaiumdar fc Mohif 120031 ILima fc iQ 
I2005D . As the largest virialized objects in the universe, 
it is thought that clusters may be less affected by un- 
known or poorly understood astrophysical mechanisms 
and it is hoped that it will be reasonably straightforward 
to compare observed catalogs with theoretical estimates. 

While many recent works have investigated the power 
of upcoming surveys and the challenges due to poorly un- 
derstood selection effects, there has been relatively less 
attention in the literature to methods of analysis. It has 
been generally assumed that the problem can be sepa- 
rated into analyzing the differential number counts as 
a function of redshift and investigating the correlations 
(either angular or three-dimensional) in the cluster cat- 
alog. Studies of galaxy surveys have traditionally used 
either the power spectrum or n point correlation func- 
tions to construct likelihood estimates, and most stud- 
ies of the galaxy cluster abundance have assumed Pois- 
son statistics. However, it has been recently emphasized 
that Poisson statistics for galaxy cl uster surveys may not 
be a particularly good assumption l)Hu fc Kravtsovl200,'^ 
lEvrard ct al. 2002). 

The importance of sample variance considerations is 
two-fold. It broadens the probability distribution of the 
counts in cells, weakening constraints on the true number 
density, and the amount of this broadening is in principle 
directly calculable from theory, adding new information 
which in principle can be used to tighten constraints. As 
will be shown below, the effect of sample variance on 
cluster counts is not negligible and treating the counts 
as a Poisson process is incorrect. 

In t his work, we will go back to the work of iWhit j 
l)1979|) in galaxy clustering and update it to the context 
of galaxy cluster studies. In particular, we will derive the 
correct likelihood function to use for binned estimates of 
the cluster number counts and we write down an explicit 
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expression for the likelihood of observing a particular 
galaxy cluster catalog as a function of theoretical models 
for the correlation structure of the density field. We will 
present some very preliminary tests using the catalogs 
from the Hubble volume simulation, demonstrating the 
regimes where these methods are robust and where more 
work clearly needs to be done. 

2. FORMALISM 

We adopt the formalism of Dodelson, Hui, & Jaffe 
(1997; hereafter DHJ97), which in turn closely follows 
the development of White (1979; hereafter W79). All 
that follows below is clearly outlined in W79 and the 
interested reader is referred to that paper for a full ex- 
position. 

The fundamental assumption is that the probability 
of observing an object at a particular location (in the 
absence of any other knowledge) is simply proportional 
to the number density of objects: 

P{Xi) = nwi{x{)dVi (1) 

where n is the number density, wi{xi) is the survey 
weight at that point (the probability of detecting an ob- 
ject if it was at that point), and dVi is an infinitesimal 
volume around the position xl. 

By similar reasoning, the probability of observing two 
objects at two points in space xi and can be expressed 
as 

PiXiX2) = n^[wiix-'i)wi{x'2) + W2ix{,x'2)]dVidV2 (2) 

where W2(a?i,a^) is the two-point correlation ftmction, 
indicating the excess probability (over random) of de- 
tecting two objects at the given positions. 
This reasoning can be extended to three objects: 

P{XiX2X3)^n^[wiix-{)wiix-'2)wiix-'3) + (3) 
W2 [Xi , Xj )wi (xl ) -I- 
W3{xi,X2,X3)]dVidV2dV3 

with W3 the three-point correlation function, and the sum 
refers to summing over all unique k. Obviously this 
can continue to even longer lists of objects. 

A given catalog will contain N objects, and the proba- 
bility of detecting those N objects at those positions can 
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be constructed from the hierarchy of correlation func- 
tions, but a very useful component of the catalog is that 
there are no objects detected at other points. A quantity 
of interest is therefore the probability of detecting N ob- 
jects at observed Xi and no other objects at the other 
positions. 

A recipe for incorporating this extra information was 
described in W79. Instead of using the functions Wi that 
were used above, one can simply replace each occurrence 
of Wi with Wi , where this latter quantity is simply 



Wi{xi,X2,...,Xi)^ (4) 

^ — — / ••• / Wi+j{xi, ...Xilj)dVi+i...dVi+j 

It is nice to have an exact likelihood formulation for a 
given catalog, but a hierarchy of infinite sums can be a 
bit unwieldy. DHJ97 noted that for galaxy surveys there 
can easily be collections on the order of 10'°^"° terms, even 
in the limit of Gaussianity and weak correlations. It will 
be shown below that galaxy cluster surveys the num- 
ber of terms is quite manageable, at least in the limit 
of the cluster distribution being Gaussian and for a rea- 
sonably sparse survey. Cluster surveys can effectively 
be considered as perturbatively different from a Poisson 
distribution, allowing the probability distributions to be 
calculated quite efficiently. 

In the Gaussian limit, it is useful to consider a few Wi: 



Wo{V) 



n 

Y 



dVidV2W2{xi,X2) (5) 



Wi{xi) ^l-n j dV2W2{xi 



X2) (6) 



W2{xi,X2) = W2{xi,X2) (7) 

where V is the survey volume and the 'Wi{xi) have all 
been set to unity for simplicity. Note that this is not 
a requirement and is a straightforward way to include 
survey non-uniformity. 

This allows an expression to be written down for ob- 
serving a collection of points as a function of a theory 
which predicts the number density and two-point corre- 
lation function. To first order in the two-point function 
(DHJ97), 



ln{P) = -nV 



(8) 



W2{Xa,Xb) 
WiiXa)Wi{Xb) 



llWiia 



In the limit of sample variance being negligible (i.e., 
in the limit where Wi approaches Wi), this approaches 
the likelihood one would expect for the number counts 
(the first term and the first term in brackets approaches 
the Poisson expression) combined independently with the 
correlation function. The second term in brackets closely 
resembles the likelihood one would obtain for the correla- 
tion function assuming that the probability of observing 
a pair at a given pair of locations is a Poisson process 
with mean given by ( 1 + W2)- When sample variance 
is not negligible, however, the probability will not be 
cleanly separable. 

Furthermore, the same formalism was used by W79 to 
write down a compact expression for the counts-in-cells 



probability: 



i-nV)^ d^ 
iV! dn" 



exp[Wo{V)] 



(9) 



where $jv(^) indicates observing N objects (and only 
N) in a volume V . 

Again, assuming weak correlations and working in the 
Gaussian limit, the counts in cells probability can be de- 
rived. Some details are derived in the appendix, but the 
result is straightforward. The probability of observing N 
objects in a volume V (to first order in the clustering) is 



N — 

X e 
iV! 



w,N-N^ N 

2^-^ + 2--^)] 
(10) 

where N is again the observed number, x = nV , and 
w = J dVidV2W2{xi,X2) is the mean correlation 

function over the survey geometry. In calculations in this 
paper we use the second order expression to generate the 
figures, but a careful analysis of terms indicates that if 
the second order term is important then there is a good 
chance that all higher order terms must be included. 

Where might this expression be expected to break 
down? In W79 it is clear that Wq starts to become a chal- 
lenging series when nVw becomes of order 0.1 or higher. 
We should therefore expect that this likelihood calcu- 
lation will be most effective in that case and that it is 
suspect outside that regime. The breakdown appears to 
arise because the assumption of Gaussianity must always 
break down: the density field cannot be negative. With 
sparse sampling the probability of encountering a neg- 
ative density is negligible and the assumption of Gaus- 
sianity is sufficient. However, at higher number densities 
the intrinsic non-Gaussianity becomes important.'^ 

Notice that this approach is straightforward to extend 
to objects that are tagged with an observable such as flux 
or richness. Going back to eqn|Hl it should be sufficient 
to replace the initial number density with a sum over 
the individual number densities in each observable bin 
(e.g., theoretical number density as a function of optical 
richness or SZ flux), and then keep in mind that each Wi 
will depend on the mass of the object or pairs of object 
under consideration. The efficacy of this technique is 
currently a topic of investigation. 

3. COUNTS IN CELLS 

We use the cluster catalogs derived from the Hub- 
ble volume simulations ( Colbcrg et al. 2000)'' to test the 
range of validity of the likelihood expression derived for 
the counts in cells. A straightforward test is to use the 
z — cluster catalogs, where evolution effects can be 
neglected. We further specialize to the tCDM simula- 
tion, where the transfer function is easy to reproduce. 
The simulation parameters are fim = 1.0, cts = 0.6 and 
F = 0.21 and a box size of 2000/i-iGpc. 

The volume is divided into subregions of cubic geome- 
try and the above expressions are tested for a variety of 
cell sizes and mass thresholds. The number of clusters 
within each cell is then used to construct the probability 
distribution of the counts in cells. Mass thresholds were 

As this work neared completion a preprint from Hu & Cohn 
appeared where they correct for this by ignoring those excursions 
of the density field and treating the remaining field as Gaussian. 
* http:/ /www. mpa-garching.mpg.de/Virgo/hubble. html 
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used to construct counts in cells for all objects above a 
threshold mass, with the mass defined as the number of 
particles (each of mass 2.2 x 10^^h~^MQ) found by the 
friends of friends algorithm with a linking parameter of 
0.2. 

The mean correlation function of the matter distribu- 
tion was derived as 



w - 



J dVidV2'W2 ix[,x'2) 

SkW{k)W*{k)P{k) 



ill) 



1 



87r3y2 

where W{k) is the window function for the cells (not 
related to the Wi above!). For cubic cells, the window 
function is simply a product of sine functions. 

The mean correlation function of the galaxy clusters 
needs to be corrected for the bias function. This is 
non-trivial; several fitting functions exist in the liter- 
ature but typical errors appear to be on the order of 
10% when applied in a regime w here they have not been 
calibrat ed iColberg et al.ll2000j) . We use the bias rela- 
tion of iSheth et al.l l)2001fl with the mean bias for the 
cluster sample above a mass threshold derived using a 
m ass function w e ighted bias factor. The mass function 
of I Jenkins erall (|2p01') was used, which has been cali- 
brated using the simulation used here. In principle, given 
a mass correlation function one can estimate the bias as a 
function of an observable cluster property and use this to 
estimate the cluster mass as a function of cluster prop- 
erties. However, current uncertainty in the theoretical 
understanding of bias as a function of mass in arbitrary 
cosmologies complicates any such undertaking. 

An example is shown in figure 1. Taking a relatively 
high threshold of 6.6 x IO^/i-^Mq and a cell size of 154 
ft.~^Mpc, the counts in cells distribution can be synthe- 
sized. We average over ten random cell offsets to smooth 
the observed distribution and compare it to the exact cal- 
culation of eqn^l In the limit of massive (and therefore 
rare and sparse) clusters the likelihood for the counts in 
cells is quite accurate. 

The assumption of a Gaussian random field will break 
down when nVw is of order 10~^ (W79). To investigate 
the nature of the breakdown, we show in figure 2 the 
first breakdown for cells of size 100 Mpc and 133 Mpc. 
At higher masses (i.e., lower number densities) the the- 
oretical result matches the simulations quite well. It is 
clear that there is noticeable deviation starting around 
nVHj ^0.5. The mass thresholds were selected to strad- 
dle the transition. At yet lower masses, the theoretical 
result breaks down fairly comprehensively, developing bi- 
modality. This is an indication that the theory is simply 
not well defined in this regime. It is impossible to set up 
a distribution of positive mass concentrations with van- 
ishing correlation functions above the two-point function. 

As the mass threshold is lowered, clusters become more 
numerous. As the cell size decreases, the mean number 
per cell will decrease but the mean correlation function 
over the cell volume will increase. For reasonable power 
spectra it turns out that the number per cell is slightly 
more important in this case for the purposes of calcu- 
lating nVw^ making it slightly advantageous to move to 
smaller cell sizes. 

To investigate the regime where this likelihood ap- 
proach can be trusted, a grid of masses and cell sizes 
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Fig. 1. — Distribution of cell counts for clusters above 6.6 X 
W^"^ Mq in Hubble Volume simulation compared to predic- 
tion from theory (broader distribution with lower peak) and cor- 
responding result from same clusters with randomized positions 
(narrower distribution with higher peak). 



was investigated. The quantity nVw was calculated and 
the lines in figure 3 indicate the locus of pairs of M and 
cell size that lead to nVw = 0.5. For a KCDM model 
with (Tg = 0.9 the likelihood approximation breaks down 
below 4 X IQi-^/i^iM©, with the tCDM breaking down 
around the same mass. However, this is a very steep func- 
tion of the amplitude of the power spectrum. For exam- 
ple, using cTg — 0.81 drops the threshold of apphcability 
to below 3 X 10^^ft-~^MQ. This is particularly relevant 
because this corresponds to the amplitude of fluctuations 
at z — 0.2 if CTg = 0.9. Any survey for galaxy clusters at 
z > 0.2 and targeting masses above a few xl0^^ft.~^AfQ 
is squarely is the regime where this likelihood calculation 
can be used. 

4. OUTLOOK 

The preliminary (and very highly idealized) tests out- 
lined above demonstrate that many cluster surveys can 
be considered to be in the regime where catalog likeli- 
hood methods will be readily applicable. Tests on counts 
in cells at z = indicate that the non-Gaussianity of the 
cluster distribution becomes important around masses 
of several times IO^^/i^^Mq. Applications of the cluster 
catalog likelihood (eqn (SJ will likely have similar prob- 
lems at lower masses and may also have difficulties at 
small separations. This is a topic of ongoing investiga- 
tion. 

The tests performed in this work used a z = simula- 
tion snapshot, rather than a lightcone selected volume. 
This was done to cleanly separate the possible complica- 
tions and there is no conceptual problem with general- 
izing this to include redshift evolution. In fact, to gen- 
eralize this one simply goes back to eqn ^3 and allows 
n, y, w to be functions of z and cosmological parameters, 
including the evolution of the relevant window functions. 
Non-cubic cell volumes that are functions of cosmologi- 
cal parameters and redshift will certainly complicate the 
analysis, but all the results presented here should be di- 
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Fig. 2. — Number of cells in simulation volume with given number of galaxy clusters above the mass threshold. The mass threshold is 
varied for two values of fixed cell size: 100 /i~^Mpc (left) and 133 h~^Mpc (right). In the left panel, the values of nViB range from 0.44 
to 0.68 (left to right in the figure) and in the right panel they range from 0.34 to 0.6 
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Mass (h-i M.) 

Fig. 3. — Combinations of mass and cell size that lead to 
nVm = 0.5 assuming Qj^j = 0.3, Q/^ = 0.7, and power spectrum 
shape parameter F = 0.21 for two values of erg. The dotted line 
corresponds to the tCDM model that was used in the Hubble Vol- 
ume simulations. The parameter space to the right of each line is 
in the regime of being easy to implement the likelihood formulation 
in this paper. 



rectly applicable. 

There are several directions for more work in this area. 
Tests on existing cluster catalogs will yield insight into 
exactly how hard it is to implement these ideas in prac- 
tice. One of the lessons learned is that this method works 
best on sparse catalogs (i.e., massive clusters), so in ex- 
isting data it would be most easily tested on the most 
massive X-ray cluster surveys. It is unfortunate that the 
method appears to break down in a regime where there 
are high quality local cluster surveys. 

On the theoretical side, an obvious next step is to apply 



these methods using clusters that have been additionally 
tagged with an observable such as X-ray or SZ flux or 
optical richness. Within this formalism such a modiflca- 
tion is entirely straightforward, as all that is required are 
number densities and correlation strengths. Incorporat- 
ing scatter in something like a mass-observable relation 
would be in principle straightforward, as would allowing 
for non-trivial (non power-law) evolution in such rela- 
tions. However, uncertainty in the bias relations as a 
function of cosmology and mass must be improved to 
realize the full potential of these methods. 

Given that the breakdown in this formalism appears 
to coincide with the onset of non-linearity, it may be 
possible to extend these results to lower masses with 
some modelling of the effects of weak non-linearity. 
The excellent match between theory and simulations at 
higher masses demonstrates that non-linearity and non- 
Gaussianity are not important in the high-mass regime. 

5. SUMMARY AND CONCLUSIONS 

In this paper we outlined a likelihood approach to 
galaxy cluster surveys that is uniquely applicable to the 
case of massive galaxy clusters. Essentially, for small 
enough sample and/or weak enough clustering strengths 
one can treat cluster catalogs as being expected to be 
drawn from a distribution that can be thought of as a 
Poisson distribution with some perturbations. By virtue 
of being only perturbatively different from a Poisson pro- 
cess, cluster catalogs lend themselves well to efficient like- 
lihood estimation techniques that naturally include both 
shot noise and sample variance. This does not require 
an arbitrary conceptual breaking of the problem into 
"counts" and "clustering" but instead simply addresses 
the question of how likely it is to obtain a set of points 
given an assumed underlying correlation structure. 

For galaxies, the number of terms involved in the es- 
timate is prohibitive, and it is also important to under- 
stand the higher order correlations. That does not ap- 
pear to be the case for massive galaxy clusters. The 
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Gaussian weak correlation limit does a very good job 
of matching the counts in cells distribution for massive 
clusters in the Hubble volume simulation. As the mass 
threshold is lower, and the number density therefore in- 
creases, some disagreement arises, likely due to the num- 
ber density becoming high enough that the distribution 
becomes sensitive to non-Gaussianity. 

It should be emphasized that it is clearly not an ad- 
vantage to have fewer clusters. One could always thin 
the number density of a sample until the catalog is in 
a regime where a catalog likelihood can be calculated. 
This would obviously be throwing away a huge amount 
of information and it seems unlikely that the gain in sim- 
plicity would be worth the loss of information. 

These techniques are straightforward to generalize to 
explicitly include inhomogeneous survey selection and 
differential number densities and clustering strengths as 
a function of galaxy cluster properties. 



With galaxy cluster surveys now being used for cosmo- 
logical parameter estimation, it is essential to include the 
effects of sample variance to obtain accurate parameter 
estimates, so obtaining accurate probability distributions 
for counts in cells should be a high priority. The tech- 
niques here are a step in that direction, showing that 
such accurate distributions are easily calculated for mas- 
sive galaxy clusters. 

As this article neared completion, an article appeared 
by Hu & Cohn discussing very similar ideas from a 
slightly different perspective. The results here largely 
agree with theirs in the aspects where they overlap. 

I have had useful discussions with Wayne Hu and have 
benefitted from an NSERC Discovery grant, the Canada 
Research Chairs program, and the Canadian Institute for 
Advanced Research. 



APPENDIX 

DETAILS OF THE COUNTS IN CELLS LIKELIHOOD 

The calculation of the likelihood for counts in cells can be derived from equation|21with straightforward algebra and 
accounting in the limit of weak clustering. Some of the details are spelled out below. 



P[^n{V)] = ^^^—exp[Wo{V)] (Al) 

Define x = nV and using Wq = —x + x^w^l'^^ the first step is to expand the exponential term as a Taylor series. 
For convenience, define a = wil1 and the counts in cells probability distribution can be written 

/ NAT JJV / -, \ 

Expanding the part inside the sum with the binomial expansion gives 

P[.,„,V,1 = <9|--i^ ^^_^,-.„.«<,> ,A3) 

\i=0 fe=0 ^ ' / 

The effect of the derivatives is that all terms in the sums with i + k < N are set to zero and all other terms are 
muhiplied by x~^{i + ky./{i + k - Ny.: 

^2=0 k=max{0,N~i) \ / \ I j 

Up to this point, the only assumption implicit in this analysis is that the correlation functions beyond the two- 
point function are negligible. To see the structure in the probability distribution, we now assume that a is small and 
investigate the relevant terms in the above equation. In particular, for small a the only important terms will be terms 
with small k. This immediately identifies the important terms as being the ones with i close to iV. 

For the terms with A; = 0, we are left with the sum 

which reproduces the Poisson probability in the case of zero correlations. 
To first order in a, one can move to one lower value of i: 

k = l: a Y.^ + = -«--^(^' -N-2XN + x^) (A6) 

and to second order in a we can continue down: 

2{i + 2 - N)\{i - 2y} ' 2 

Ax^N + IIN^ - 6N^ -6N - ixN^ + x^ + N^) (A7) 



- E ,U ^ O -^ U. - .V (-)"^'" - ^e-^-,xN + 12.A.- ,x^N^ - 6x^N ■ 
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with subsequent orders obtained similarly. These sums are straightforward to implement using software such as Maple 
but can become tedious. In all that was implemented in this work the series was truncated at second order. 

The sums have an interesting property that each term seems to be of order {ax^)^ , since A'' is of the same order as 
X. At each level the terms are actually of lower order in x due to cancellations in the sums. For example, for x — N 
the expression in brackets in equation IA7I reduces to 37V^ — QN , making this term roughly a^x"^. Nonetheless, it is 
clear that simply having small a is not sufhcient. For this formalism to be useful it is required that ax be small, at 
least, and it is preferable that ax"^ be small to be sure that all higher order terms are negligible. It is apparent from 
the expression for Wq that bad things are happening as ax is becoming large, since it is becoming dubious that the 
void probability function is bounded to be less than one (i.e., Wq must be negative). This is likely a breakdown of the 
physical picture of Poisson sampling of a Gaussian random field, rather than a breakdown of any approximations. It 
is impossible to Poisson sample with a negative mean, which a Gaussian random field allows. The requirement that 
the density is non-negative necessitates some amount of non-Gaussianity. With enough samples, this non-Gaussianity 
must manifest itself in the observed distribution. 
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